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The solutions to the heat diffusion equation for spherical absorbing centers in laser materials are 
presented. The ratio of the volume specific heats for the absorbing center and the laser host is found 
to determine three regions of behavior. Series expansions for small times and for very large times also 
are given in each of the three regions. Rapidly converging representations for the complementary ('nor 
function of complex arguments are developed in order to evaluate numerically the region for which 
the volume specific heat ratio is greater than three fourths. 
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1 . Introduction 

In a previous paper [1], a model to describe the 
behavior of absorbing centers in laser materials has 
been developed. Computing the thermal stresses pre- 
dicted by this model requires solutions to the heat 
diffusion equation for the temperature. The mathe- 
matical details for solving the heat diffusion equation 
are omitted in ref. 1 and are reported in the present 
paper. In particular, solutions to the heat diffusion 
equation for spherical absorbing centers in laser ma- 
terials, series expansions of these solutions at the 
center-host interface when times are sufficiently small 
and large, and a proof that in the limit of small inclu- 
sions the temperature of the center-host interface is 
independent of the bulk thermal properties of the 
inclusion are presented. These solutions may be 
applied also to other diffusion problems such as oxygen 
diffusion in Ti0 2 , AI2O.3, ZnO, and other oxides [2]. 

When certain physical conditions exist, some solu- 
tions to diffusion equations contain error functions for 
complex arguments. For example, when the volume 
specific heat of a spherical absorbing center is greater 
than three fourths times the volume specific heat of the 
host, the analytic expressions for the temperature in 
the host contain error functions for complex arguments. 
This is the case for platinum inclusions in some neo- 
dymium doped laser glasses. The oxygen diffusion 
problem is a second example. When the volume diffu- 
sion coefficient is sufficiently greater than the surface 
exchange coefficient, the analytic expressions for the 
oxygen concentration profile also contain error func- 
tions for complex arguments. This may be the case for 
oxygen diffusion in zinc oxide [3J. Representations of 
the error function for complex arguments which yield 
efficient computer programs are not readily available 
in the literature and therefore are derived in the 
appendix. 



2. Mathematical Description of the Model 

The heat diffusion equation and the boundary 
conditions are summarized in this section. 

A spherical absorbing center of radius r (cm) is 
imbedded in a host of infinite extent. The absorbing 
center and host are initially at a uniform temperature 
7V An incoming spherical wave of light (laser beam) 
with an energy flux (7rr^EiJ^7rr 2 r) (W/cm 2 ) falls at 
time £ = upon the absorbing sphere. The quantity 
El is the energy density (J/cm 2 ). The quantity r=|r| is 
the distance from the origin of the system, which is 
the center of the sphere. A square-wave light pulse of 
width t(s) is assumed. The sphere absorbs an equiva- 
lent energy flux //(f) (W/cm 2 ) uniformly over its entire 
surface. The energy flux for a square-wave laser pulse 
is 

f 0, when t < and t > r, 
H(t)=\ (1) 

A(\, T„)(EJ4t\ when ^ t ^ r, 



where A(k, T a ) is the absorptance at absolute tempera- 
ture T a for incident radiation having wavelength \. 
The problem is to compute the temperature 0(r, t) of 
the absorbing center-host system as a function of 
time t and distance r= |r| . 

The heat diffusion equation describes the time and 
space dependence of the temperature for large times 
t and large distances r; namely, when r < r it is 



(a0(r, t)/dt)=a*.V 2 e(r, t), 
and when r > r () it is 

(a6(r, t)/d$)=atV*e(r 9 t). 



(2) 



(3) 
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The Laplacian operator is denoted by V 2 and the 
thermal diffusivity is a 2 = (K/pC) (cm 2 /s), where K is 
the thermal conductivity (W/cm °C), p is the density 
(g/cm 3 ), and C is the heat capacity at constant volume 
(J/g °C). The subscript c refers to the absorbing 
center and the subscript h refers to the host. The 
diffusion equations (2) and (3) are valid only when 
local thermodynamic equilibrium exists. They become 
suspect for times comparable to the phonon-phonon 
collision time ( ~ 10~ 12 s) and for distances comparable 
to atomic dimensions (~ 10~ 8 cm). The diffusion equa- 
tions require a statement of the boundary conditions 
before solutions are uniquely defined. The temperature 
has the form, 

f To, t^O 

8(r,0= (4) 

[ T + T(r 9 t), t>0 9 

for all values of r. The function T(r, t) is the tempera- 
ture for t > relative to the ambient temperature 7V 

The boundary conditions are time dependent. The 
temperature must be finite everywhere. 



3. Laplace Transform Solutions 

Equations (8) and (9) and the initial conditions and 
boundary conditions (4), (5), (6), and (10) are solved by 
taking their Laplace transforms. The Laplace trans- 
form of the temperature F(r, t) is denoted by 



U(r,s) 



i 



exp (-st)T(r, t). 



The Laplace operator J? is introduced and is defined 
bjJ?{T(r, t)} = U(r,s). 
The Bromwich integral, 

"I C c+id 

F(r,0=Lim— exp (st)U(r,s)ds, (11) 

d^ x iTTl Jc-id 

expresses the temperature in terms of the Laplace 
transform. The real quantity c is chosen to be suffi- 
ciently large so that the integral, 



/; 



exp (—st) \T(r, t) \dt, 



(12) 



r c (r=0, t) = finite and Th(r=°°, t) = finite. 



(5) 



The continuity of the temperature across the center- 
host interface is 



T c (r = r r, t) = Th(r=r r, t) , 



(6) 



where r is the radial unit vector. The conservation of 
heat flow at the interface gives 



exists. The inverse Laplace operator J? -1 is defined by 
o§? -1 {£/(r, s)} = T(r, t), and is a convenient way to 
represent the integral in eq (11). The Laplace transform 
of the energy flux H(t) is denoted by h (s). Taking the 
Laplace transform of eq (1) yields that 



h(s) = (QI4rs){l-exp(-sT)}, 



whe 



(13) 



-K c VT c (r=r Q r, t) -r + H(t) Q = A(\,T a )(E,J4) 

= -K h <7T h (r=r r,t) • r. (7) 

The Laplace transforms of eqs (5), (6), (8), (9), and 
Equations (1) through (7) represent the mathematical (10) are respectively, 
description of the model for absorbing centers in laser 



materials. The physics from which these equations 
arise is stated in ref. [1]. It is most convenient to use 
the spherical coordinates r, 6, and <p because the heat 
flux H(t) is independent of 6 and (p. The problem then 
becomes a spherically symmetric one. Equations (2), 
(3), and (7) now have their respective forms, 



ar c (r, t) 



d 2 T c (r, t) , 2dT c (r, t 



dr* 



+ -■ 
r 



dr 
for 



r < r , 



dT h (r,t) 



and 



r, t) _ ( d*T h (r 9 t) 2 dT h (r,t) \ 

t a/, \ dr* r dr /' 



K c dTc{r ~ r(ut K H(t)=-K h 
at 



dr 
for r > r , 

dT h (r=r , t) 
dr 



(8) 



(9) 



(10) 



t/ c (0)= finite and £/„(») = finite, (14) 

U c (zco)=U h (z M ), (15) 



cPU c 



dU c 



dz f +2z ^-^ + 



7 , (r,0) z 2 



d*U h 



dU h 



'* j 2~ +2z »~i — 4 U » + 

azi dzh 



To(r,0)zi 



— K c -: h h(s) 

a c dzc 



K h 



= for r < r , 
(16) 

= for r > r , 

(17) 

s l > 2 dU h (z h = Zho) 
a h dzh 

(18) 
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where 



Zc= (rs ll2 la c ),z h = (rs 1 l 2 la,,),z C ()= (roS ll2 /a c ) 

z h0 = (r {) s^la h ),U c (r,s) = U c (z c ), 

and U h (r, s) = U h (z h ). 

General solutions to the Laplace transformed diffu- 
sion eqs (16) and (17) are 

z z s 

Let us for convenience choose the zero of the tem- 
perature scale such that To = 0. The condition, 
U c (0) = finite, requires that B> € = and gives 



U c (zc) =Bic{sinh {z c )lz c }. 



(20) 



The condition, Ut, (°°) = finite, requires that Bit, -+- B%h =0 
and gives 

U f ,(z„) = B 2h {exp (-Zh)lz h }. (21) 

The continuity condition (15) relates the constants 
B\ c and B 2 i,; namely, 

B u = B 2 h{z c o exp (— Zho)lzho sinh (z c o)}. (22) 

Finally, the conservation of heat flow through the inter- 
face, eq (18), and the condition (22) permit us to express 
the constant B>t, in terms of r<>, K r , a c , Kt,, and a*; 
namely, 

# 2 » - h(s)(z h0 ls 1 ' 2 ) exp (zho) X [(ftk/o*)(l + Zno 1 ) 

- (Kd a c )(z^-coth (zco)]" 1 . (23) 

Substituting eq (23) into eq (21) and evaluating the 
Bromwich contour integral to obtain the temperature 
of the host Tt, (r, t) would be a most ambitious task. 
Two different arguments exist by which eq (21) reduces 
to a more tractable form. The first argument treats 
eq (23) in the limit that z c o < 1. Then the coth (z c o) in 
eq (23) is approximated by 

COth (Zco) ~ Zco 1 + (z C ()/3) + W c o) + • • - , 

and eq (21) thereby is approximated by the expression 
U hx (z h ) = Hm Uh(zh) = B„(oo, <>){exp (- zh)Izh}, (24) 

Zctr+0 



where 



fi/,Ks) = 



h(s)zho exp (zho) 



Jlffa^ + fc+M^ + M' 
M= (pcC c ro/3), />. =b{\T (l-fl)i/*}, 
6 - (K h l2Ma h ) , and /? = (4p c Cc/3p ft C/, ) . 



The second argument involves simplifying boundary 
condition (10). Observe that when £ r ^> /£/, the tempera- 
ture gradients in the absorbing center are much less 
than those for the host. In the limit K r $> Kn boundary 
condition (10) becomes 



M JJj^= Kh >J^A +H(th 



dt 



dr 



(26) 



where r^ = Lim (r ±e). The temperature of the sphere 

e-*0 

becomes spatially uniform; that is, T c (r, t)=T c = Th{r {) , t) 
whenever its thermal conductivity is much larger 
than that of the host. The heat flux into the sphere is 
given by, 

{ (47rrg/3)pcCc(dr c /^)/47rrg} = M(dTJdt). 

Taking the Laplace transform of the conservation of 
heat flux, eq (26), which is valid when K c >Kh, the 
author obtains the same expression for the Laplace 
transform Uh as that given in eq (24). Hence, approxi- 
mation (24) obtains either when |r 5 1/2 /a r | < 1 or when 
K c > Kh. The latter inequality imposes no restriction 
on ro and s or equivalently on r {) and t but the former 
inequality does impose such restrictions. 

The volume specific heat ratio R = (4p c C c /3p/iC/i) 
determines three regions of behavior. Namely, when 
/? < 1, the roots b± are real and are not equal; when 
J? = l, the roots are real and equal; and when R > 1, 
the roots are complex and conjugate to one another. 
The latter is the region for which the arguments of 
the complementary error functions are complex. 
In the next Section, analytic expressions for the 
temperature Tt, (r, t) in each of these three regions are 
presented. 

4. Results and Conclusions 

In this section, closed form expressions for Tt,(r, t) 
are given for the regions R < 1, /? = 1, and R > 1. 
Also, series expansions for large and small values of 
bt 112 are derived. 

In the limit that z c0 ^ 1 or in the limit that K ( > Kh, 
the Laplace transform of the temperature in the host 
becomes 

U hx (z h ) = [l-exp(-ST)]D h (s) (27) 



where eqs (13) and (25) require that 

(Q/4ts)z;i exp (zhp—Zh) 



D h (s) 



2/,M(5 1 ' 2 + 6 + )(5 1 / 2 + 6-) 



(28) 



The translation property of Laplace transforms, (eq 
(29.2.15) of ref. [4]) then permits one to express the 
temperature in the host in terms of a Bromwich inte- 
gral, Bi(r, t), and the unit step function, u(t); namely, 
[5] 



T h (r, t)=B,{r, t)-u(t-T)B,(r, t-r). (29) 



(25) where 
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£ 7 (r, t)=Lim 



\ Cc + id 

tri Jc-id 



exp (st)D h (s)ds, (30) 



and the unit step function u{t) is denned by 

0, t<0 
u(t)= i, t = 

1, t>0. 



(31) 



Subsections 4.1, 4.2, and 4.3 of this section contain 
respectively the evaluation of the Bromwich integral 
(30) for the cases of R < 1, R = 1, and R > 1. 

4.1. R< 1 

When the volume specific heat ratio R is less than 
one, then the roots b + and 6- are real and are not equal. 
The Bromwich integral (30) is evaluated more readily 
by expressing D h (s) in terms of partial fractions; 



where 



D h (s)=P<(r){d + (s)-d-(s)}, 
P<(r) = {QnSrMbr(l-Ryi 2 }, 



(32) 



(33) 



and 



d>(s) = [exp{(r {) -r)s^la f) }ls(s^ + b,)] (34) 



The transform (34) is valid only when r > r . The 
Bromwich integral (30) is evaluated by applying eq 
(29.3.89) of ref. [4] to eqs (32) and (34). This evaluation 
is correct only when y = {(r—r )la h } ^ 0. It is found 
found that according to ref. [4], 

T h (r,t;R<l) = P<(r)[b-i{F(y,b + ,t) 

-u(t-T)F(y,b + ,t-T)}-bz l {F(y,b-,t) 

-u(t-T)F(y,b-,t-r)}] (35) 

where y= {(/ ' — r )/a/j} ^ 0. The dimensionless func- 
tion F(y, 6, t) contains the complementary error 
function; namely, 



F(y,b,t) 



where 



exp (yb+bH) eric(bt 1 ' 2 + (y/2t 1/2 )) 

+ erfc(y/2r 1 / 2 ) 



erfc(x) = l %l | exp (-t 2 )dt. 

T7 1 ' 2 Jo 

4.2. R = l 



(36) 



(37) 



and the Laplace transform factor Di,(s), eq (28), 
becomes 



D h (s)=Po(r)do(s), 
rVr) = (Oo/4rM& 2 r), 
d (s)=b 2 {exp (-y»«*>/a(*i/* 
Because the identity, 

»)) 

exp (— ys 1/2 ) ] _ d 



where 
and 



, d I exp (-ys 1/2 ) ' 
1 <M 5(5" 2 - 



b) 2 }. 



(38) 
(39) 
(40) 



= =Sf- 



•}-; 



{b- l F(y,b,t)}, (41) 



s(a" 2 + &) 2 J db 
exists, the temperature in the host has the form 
T h (r,t;R = l)<=Po(r){G(y, b, t) 

-u(t-T)G(y, b, t-r)}, 



where 
G(y,b,t) 



(42) 



b*-^{b-*F(y,b,t)} 



(yb + 2b 2 t+l) exp (yb + b 2 t) erfc (bt 112 
+ (y/2f 1 ' 2 )) + erfc (y/2f" 2 ) 

-2(6t 1 / 2 /77 1 / 2 ) exp (-y 2 /4t). 

4.3. R>1 



(43) 



When the volume specific heat ratio is greater 
than one, the roots are complex and conjugate to 
one another; that is, b ± — v z hifi where v = b and 
fi = b(R - 1) 1/2 for R > 1. Because b- = 6 J, the Laplace 
transform factor Dh (s) [eq (28)] becomes 



^here 



D h (s)=P > (r)d > (s) 



P>(r) = (Qrol4rMbrfjL), 



(44) 



(45) 



and 



d>(s) 



_ exp ( — ys ll>2 ) 
lis 



1 



1 



|l/2 + l/- 



+ v—in s ll2 -\-v-\-ifx\ 



(46) 



Evaluating the contour integral which results after 
eq (44) is substituted into eq (30) by the calculus of 
residues is a most difficult procedure. Because the 
integrand in eq (37), exp (— £ 2 ), is analytic everywhere 
When the volume specific heat ratio R equals one, in the finite complex plane, the function, b~ l F(y, b,t), 
the roots are real and equal; that is, may be analytically continued to complex values of 

b. The temperature then is computed by taking 
b + = b- = b= (Khl2Mcih) the imaginary part of the analytic continuation of 
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6 _1 F(y, /;, t) for complex values of b; namely, 

T h {r,t;R>l)=P > (r){H(y,v,9i,t) 

-u(t-r)H(y,v,iJL,t-T)}, (47) 

where 

H(y, v,fi,t) = b Im {(*> — i/A) _1 F(y, v — ijx, t)} 



The analytic expression for the function H(y, v, p,, t) 
is rather lengthy. It contains factors related to the 
real and the imaginary parts of the complementary 
error function for complex arguments. The error 
function for complex arguments and procedures for 
its numerical evaluation are given in Appendix A. Let 
us introduce the following notation: = ix(y + 2vt) , 
p= vt xl2 + (y/2f 1/2 ) and a> = pf 1/2 . It is shown in 
Appendix A that the real and the imaginary parts of the 
complementary error function for complex arguments 
contain the following integrals: 



•»/: 



sdaw (p, co) = — : — exp (— o> 2 ) J exp (£ 2 ) sin (2p£)df 

and 

2 f w 

cdaw (p, (o) = —. —- exp (— co 2 ) exp (£ 2 ) cos (2p£)d£. 

TT 1 '" J<> 



The expression for H(y, v, p, t) becomes, in terms of 
this notation: 

#(y, ^,/x,0 = {W(^ 2 + i^ 2 )}[^erfc (y/2^/ 2 ) 

— exp {yv + v 2 t — fjL 1 t){{—v sin 0+ p cos 0) 

X (erfc (p) —exp (a> 2 — p 2 ) sdaw (p, w) ) 

+ (i^cos 0-f psin 0)exp (co 2 — p 2 ) cdaw (p, co)}]. (48) 



At first one may question the oscillatory factors, 
cosine of and sine of 0, which appear in eqs (47) 
and (48). One does not usually expect such behavior 
for solutions to a diffusion equation which has no 
oscillatory boundary conditions or source terms. How- 
ever, the oscillatory factors in eq (48) are over damp- 
ened so that the spatial gradient is always less than 
zero for all r ^ r [6]. Similar oscillatory factors arise 
in other diffusion problems. Theoretical solutions for 
fluid turbulence and for the characteristic lines pro- 
duced by a moving disc in a uniformly rotating, non- 
viscous, and incompressible liquid also may contain 
oscillatory factors [7J. The oscillatory factors in eq (48) 
exhibit a wavelength k osc = (27ra/,/p,) , a frequency 
(1/Tqsc) = (ixv/n), and a velocity v osc = — 2vahlr. 
The oscillatory speed, |v osc | = 2va\ x is less than the 
speed of sound in the host for all the situations dis- 
cussed in ref. [1]. If this were not the case then the 
quasi-static linear thermoelastic theory would be 
suspect. 



4.4. Expansions for bf V2< ^l 

In order to avoid overflow complications with the 
computer it is necessary to use series expansions for 
the functions when they are evaluated at the center- 
host interface, y = 0. These expansions readily show 
that as bt 112 approaches the surface temperature 
Th(ro, t) depends only on the properties of the ab- 
sorbing center. This agrees with the intuitive state- 
ment that for sufficiently short times very little heat 
enters the host. The expansions for the F 9 G, and H 
functions when bt l/2 <$ 1 are listed now: 

F(0, 6, t) ~ (2l7T i ' 2 )bt'l 2 - (^/ 2 ) 2 + (4/3w 1 ' 2 )(fo 1 '*) a 

+ ^(^ 1/2 ) 4 + . . ., (49) 

G(0, 6, t) - (bt^ 2 ) 2 - (8/3tt 1 / 2 )(^ 1 / 2 )^ 

+ e(bt l i 2 y+ . . ., (50) 

and 



//(0, v, //,, t) ~ b(Mt[l - (8/3tt 1 / 2 )^ 1 / 2 



•I 



(51) 



where v = b. 

Finally, substituting eqs (49), (50), and (51) respec- 
tively into eqs (35), (42), and (47) yields when bt 1/2 <tl 
and t ^ r, 



T h (r ,t) ~«fc/4rJ#). 



(52) 



The quantity M depends only on the properties of the 
absorbing center. 



4.5. Asymptotic Expansions for large bt 12 >] 

In this subsection, the asymptotic expansions for 
F. and G when bt l/2 >\ are listed below. Such ex- 
pansions for H are not as straightforward as those for 
F and G. In fact, the derivation of a representation 
for H which converges rapidly when vt xl2 > 1 is most 
essential in the development of a computer program 
to evaluate T fl (r, t) for the case in which R>\. 
Representations for H are developed in appendix A. 
These expansions show that as r () approaches the 
surface temperature Th(r , t) depends only on the 
absorptance of the interface and the properties of 
the host and it becomes independent of the other 
properties of the absorbing center. The asymptotic 
expansions for the functions F and G when bt xl2 > 1 
are respectively, 

F(0, 6, t) ~ 1- (Hir^bt 1 ' 2 ) + . . ., (53) 



G(0, 6, t) ~ 1- (2l7T'l 2 bt" 2 ) +L . .. (54) 

Substituting eq (53) into eq (35) gives expressions for 
the interface temperature when r is very small; 
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namely, if t > r and br 1 ' 2 > 1, then 

T h (r ,t)~(Qr 2 l87T^a h K h t*lz), 
and if t = r aaidbr 112 > 1, then 

T h (ro 9 T)-(Qr l^rK h ). 



(55) 



(56) 



Observe that eqs (55) and (56) do not contain any of 
the bulk thermal properties of the absorbing center. 
The expressions for Th (r, t) given above form the 
mathematical framework upon which the computer 
program used to generate the tables of ref. [1] is 
developed. In addition, several first and second order 
derivatives of the functions F, G, and H with respect 
to time and space are needed in order to evaluate the 
temperature for large times and in order to check 
numerically that the temperature gradient is always 



negative and that the rate of expansion of the sphere 
never exceeds the speed of sound in the host. Because 
closed form expressions for the derivatives are ob- 
tained by conceptionally straightforward methods and 
because these expressions are very lengthy in the case 
of the H function, the author chooses not to list them 
here. Once the temperature has been calculated then 
the optical path length, the stress components, and the 
stress birefringence near the absorbing center are 
readily found. These problems have been discussed in 
detail in ref. [1]. 
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5. Appendix A. Complementary Error Function for Complex Arguments 

It is convenient to define 
I(x, y) = cdaw (x, y) + i sdaw (x, y) 



A computer program which evaluates the integral, 
T„(p, t)p 2 dp, for R>1, 



.1/2 



/: 



exp (— y 2 ) I exp (s 2 + 2ixs) ds , (A3) 



where 



hundreds of times must incorporate a rapidly converg- 
ing representation of the complementary error func- 
tion, if the total computation time is to be reasonable. 
Several representations of the complementary error 

function and those representations which are used in c w \ x *y) ~ ^ e i*\ x > y)s 
the studies on inclusions in laser materials [1] are 

presented in this Appendix. Computer programs for — _z_ exp (— y 2 ) e xp (s 2 ) cos (2xs)ds 
all of the representations given here have been vr 1/2 Jo 

written in Fortran [8J. j 

Let us begin with the analytic continuation of the 
complementary error function to complex arguments §daw ( } =Im , /( )} 
z = x — ty, where x and y are real, lhe complementary 
error function for real arguments z = x is defined by 9 fy 

= —r^ exp (— y 2 ) exp (s 2 ) sin (2xs)ds. 
2 fx 7T ! Jo 

erfc (x) = 1 777 exp ( — t 2 )dt. 

7T I/Z Jo 



(A4) 



(A5) 



Because the integrand exp ( — t 2 ) is analytic every- r (z) = erfc (x) 

where, the value of the integral is independent of the 
path along which the variable t follows in going from 
t = to t = z= x — iy. The following path is chosen: 



Inserting these definitions into eq (A2) yields 
: (z) = erfc (x) 
— exp (y 2 — x 2 ) {sdaw (x, y) — i cdaw (x, y)}, 



2 f x 
erfc (z) = 1 tt: exp { — t 2 )dt 

7T 112 Jo 

~^mj x ex P (~t 2 )dt. 



where z = x — iy. The definition, 



(A6) 



(Al) 



erfs (z = x — iy) = erfc (x) — exp (y 2 — x 2 ) sdaw (x, y), 
is introduced. 

Many rational approximations for the error function 
are listed in ref. [4], sections 7.1.25 to 7.1.28, and are 

Letting t = x-is in the second integral, the author readily modified to evaluate exp (* 2 ) erfc (*). One 
• . possibility uses the rational approximation section 

7.1.26 of ref. [4]. It is the representation used in ref. 
er f c ( z ) = er f c ( x ) [!] whenever z is real and z = x < 2.1; namely, 



2i f y 

H 777 exp (—x 2 ) exp (s 2 -\-2ixs)ds. 

TT l Jo 



(A2) exp (* 2 ) er ^ c (*) = 2 a >^ w + ex P ( x2 ) € (x)-> 
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(A7) 



t= (l+px)- 1 and \e(x)\ 



1.5 X10- 7 . 



cos (2rnrsla) 



(a/vr 1 / 2 ) 



1+2^ 



exp (— n 2 a 2 ) cosh (2nsa) 



(A13) 



The numerical constants p, ai, a2, #:*, cm, and ar> are 
given in ref. [4] section 7.1.26. The last term estimates 
the error and is large for large values of x. 

Two series expansions for small values of x are where a is a real positive number, < a ^ 1. The 

constant a is chosen so that 



2 x f—l\n x 2n+l 

erfc(*) = 1-^2 n! ( 2n + 1) ' 



(A8a) 



£ = 2 2 exp (-n 2 Tr 2 la z ) cos {2nns/a) « 1. 



2 °° 2 n x 2n 
erfc (a) = 1 jyj exp (-* 2 ) ]T — — 



When a = (1/2), the term E is of the order of 10~ 17 
and eq (A13) becomes an approximate expression for 



(2rc+l) exp (s 2 ); namely, 
(A8b) 



An asymptotic expansion for large values of x exists; 
namely, 



tt^ ,2 x exp (x 2 ) erfc (x) ~~ 1 



+ 1 



(— l) w l - 3 . . . (2n-l) 

(2x 2 )' 1 



(A9) 



exp (s 2 ) ~ (a/Tr 1 / 2 ) 

OC 

1 + 2 Y exp (— n 2 a 2 ) cosh (2rasa) 

+ 0(2 exp (-7r 2 /a 2 )). (A14) 



The integrations over s in the eqs (A4) and (A5) can be 
performed in a closed form when approximation (A14) 
Equations (A8) and (A9) may be continued analyti- replaces exp (s 2 ) in the integrand. 

cally to complex values z = x — iy^r exp (— i/3), where 

r= (x 2 + y 2 ) 1/2 and /3=tan~ l (ylx). The series expan- • 

sions for small \z\ = r becomes 



r ,v , 2 A (-l)"r 2 » + '{cos ((2/i+l)j3)-isin ((2n+l)j3)} 
UJ ^ 2 £ n!(2n+l) 



(A10) 



erfc (2) = 1 



grexp (~x 2 + y 2 ) * 2"r 2 "{cos ((2n + l)/3 + 2*y) - i sin ( (2n + l)p + 2xy) } 



.1/2 



X 



1-3 .. . (2/1 + 1) 



(All) 



The asymptotic expansion for large \z\ and for 
|/3 1 < (3ir/4) is 

erfc (z)~ CXP ( ~^ + y2) {cos (/3 + 2xy) + i sin G8 + 2xy)} 



+ 2 



(-l)"l-3 . . . (2/1-1) 

2 n r 2n 



{cos ((2n-l)/3 + 2*y)-isin ((2/i — l)/3 + 2*y)} (A12) 



Representations for the functions cdaw (x 9 y) and The results are 



sdaw (x, y) follow from eqs (A6), (A10), and (All), and 
(A12). 

A corollary of Poisson's formula is the basis for 
another approximation for the complementary error 
function. It is more convenient in this case to first 
approximate cdaw (*, y) and sdaw (x, y) . One may 
derive by means of Poisson's formula (7) the identity, [9] 



cdaw (x, y) jtj exp (— y 2 ) 



a sin (2xy) a ^ exp (— n 2 a 2 ) 



„l/2 2j 



exp (s 2 ) 



1 + 2 £ exp (- n 2 77 2 /a 2 ) 



2x ' tt" 2 £ t (n 2 a 2 +x 2 ) 
(na sinh (2nay) cos(2;ty) 
+ x cosh (2nay) sin (2xy)) 



(A15) 
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sdaw (x, y) 



.1/2 



exp (-y 2 ) 



a (1 — cos {2xy) ) a ^ exp (— n 2 a 2 ) 
r l / 2 ~ 2x ir^i (n 2 a 2 + x 2 ) 



(x + na sinh (2nay) sin (2xy) 

— x cosh (2nay) cos (2xy)) [• (A16) 



Among the above approximations which are valid 
for large | z | approximations (A15) and (A16) converge 
most rapidly. However, they are still not as efficient 
as the algorithm for the real and imaginary parts of 
the function, 



JF(z)=exp (— x 2 ) erfc {—ix), 



(A17) 



which is given in ref. [10]. The algorithm is based upon 
the continued fraction representation, 



exp (z 2 ) erfc (z) = 



J 1 (1/2) 1 (3/2) 2 

77 1 /' 2 Z+ Z+ Z+ Z + Z + ' ' " 

(A18) 



and efficiently separates the real and imaginary parts 
ofeq(A18). 

The complex function I(x, y), eq (A3), can be ex- 
pressed in terms of the function W{z); 

I(x, y) = i exp (x 2 — y 2 ){ exp (— x 2 )W{ix) 

-exp(-(x-iy) 2 )W(ix + y)}. (A19) 

Introducing the notation W(ix) = W\ + iW-> and 
W(ix + y) = W 3 + iW 4 , where W u W^ W^ and W A 
are real, produces the result 



cdaw (*, y) 



and 



W 2 exp (- y 2 ) + { W 4 cos (2*y) 
+ r 3 sin (2xy)}, 



sdaw (x, y) = SFi exp (— y 2 ) + {JF 4 sin (2xy) 

— W% cos (2xy)}. 



(A20) 



(A21) 



A Fortran program based upon the algorithm of ref. 
[8] to evaluate the real and the imaginary parts of the 



function W [i.e., Re {W{z)} and Im {W(z)}] has 
been written. 

Computer programs have been written to evaluate 
all of the above representations for the complementary 
error function. The prescriptions given below are more 
efficient than the many others examined during the 
computations involved in ref. [1]. First consider the 
case for lm{z} = 0; that is, z = x and y=0. The ra- 
tional approximation (A7) is used when 0<|oc|^2.1 
and the continued fraction representation eq (A18) 
is used when 2.1 <|x|< °°. This combination is effi- 
cient for study of inclusions in laser materials when 
the volume specific ratio satisfies the inequality/? ^ 1. 
Next consider the case when z = x — iy. The series ap- 
proximation (A10) is used whenever both 0< \x \ ^ 0.001 
and 0<|y|^0.001 and the algorithm of ref. [8], is 
used for all other values in the complex plane. This 
combination is efficient for the study of inclusions in 
laser materials when the volume specific heat ratio 
satisfies the inequality/? > 1. 
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